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Abstract 


A  reanalysis  method  to  analyze  the  strength  of  a  structure  which 
has  encountered  a  laser  strike  is  developed.  The  method  accounts  for 
the  following  types  of  laser  induced  damage:  1)  loss  of  structure  due 
to  melting;  2)  change  of  material  properties  due  to  temperature  changes; 
3)  addition  of  load  due  to  thermal  stress.  The  program  uses  heat 
balance  calculations  over  successive  finite  time  increments  on  an  array 
of  finite  elements  bisecting  the  laser  beam  spot  to  determine  the  tem¬ 
perature  distribution.  These  results  are  then  converted  to  structural 
stiffness  parameters  and  the  structural  analysis  is  performed  using  a 
finite  element  based  reanalysis  method.  The  reanalysis  method  predicts 
the  damage  effects  from  the  initial  undamaged  solution.  The  program 
was  found  to  give  good  results  consistent  with  results  obtained  from 
a  separate  analysis  for  each  damage  condition,  but  with  less  computer 
time  and  manhours. 


AFIT/GAE/AA/83M-1 


REANALYSIS  METHODS  FOR  STRUCTURES 
WITH  LASER  INDUCED  DAMAGE 

I  Ini  faction 

With  the  current  thrust  towards  laser  technology  in  the  development 
of  weapons,  methods  must  be  developed  for  analyzing  the  damage  due  to 
lasers.  The  total  problem  is  a  merging  of  two  separate  technologies. 

The  first  is  the  characterization  of  laser  damage.  The  thermal  problem 
occurring  in  the  area  local  to  the  beam  must  be  related  to  such  concepts 
as  stiffness  which  are  required  for  a  structural  analysis.  The  second 
is  the  area  of  structural  analysis.  A  cost  effective  method  for  per¬ 
forming  the  structural  analyses  necessary  for  analyzing  each  damage  case 
must  be  developed. 

Laser  damage  studies  have  been  conducted  to  calculate  the  tempera¬ 
ture  distributions  and  "melt-through  time"  for  specific  conditions 
(Ref  1-4).  The  primary  interest  of  these  studies  was  to  determine  the 
material  degradation  in  the  immediate  vicinity  of  the  beam  spot.  The 
object  of  the  present  study  is  to  relate  these  local  damage  mechanisms 
to  a  set  of  global  variables  which  can  be  used  in  a  structural  analysis. 

Depending  on  which  aircraft  structural  systems  are  affected  and  the 
degree  of  the  damage,  the  effect  on  the  performance  of  a  damaged  air¬ 
craft  can  vary  from  minor  changes  to  the  loss  of  the  aircraft.  Since 
military  aircraft  are  expected  to  encounter  damage,  it  Is  Important  that 
they  retain  some  level  of  structural  Integrity.  This  must  be  provided 
for  during  the  design  phase.  A  simplistic  approach  might  be  to  Increase 
the  sizes  of  all  members  and  components  which  make  up  the  structure. 
However,  the  constant  quest  to  reduce  weight  and  cost  while  increasing 
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performance  through  optimization  is  at  odds  with  such  an  approach, 
unless  limited  to  a  minimum  increase  of  critical  members.  The  problem 
lies  in  determining  the  location  of  the  critical  members,  for  It  Is 
impossible  to  predict  In  advance  the  degree  and  location  of  the  damage. 

Thus  an  analysis  must  be  conducted,  investigating  a  large  number  of 
possible  levels  and  locations  of  damage  to  identify  the  areas  needing 
modification.  Repeating  a  complete  structural  analysis  for  each 
hypothetical  damage  case  would  be  extremely  costly. 

Recent  studies  have  developed  Iterative  methods  to  reduce  the  cost 
of  the  large  numbers  of  analyses  (Ref  5-9).  Such  techniques  have  also 
been  used  to  analyze  the  effects  of  conventional  weapons  (Ref  9).  The 
object  of  the  present  study  is  to  modify  and  extend  the  reanalysis  itera¬ 
tion  technique  developed  by  Venkayya  (Ref  5)  to  allow  for  efficient 
analysis  of  a  structure  subject  to  many  different  damage  conditions,  and 
to  apply  this  methodology  to  the  analysis  of  laser  damage. 

Combining  the  objectives  of  the  two  merging  technologies,  the  total 
program  objective  Is  to  develop  a  method  by  which  a  structural  designer 
can  determine  critical  members  of  the  structure  needing  strengthening  or 
redundancy  to  survive  laser  strikes  at  relatively  low  cost.  This  paper 
follows  the  logical  development  of  the  program:  Section  II  discusses 
the  structural  analysis  method  to  be  used  and  develops  the  iterative  re¬ 
analysis  technique  to  be  applied  to  the  damaged  structure.  Section  III 
presents  the  solution  to  the  heat  conduction  problem  resulting  from  the 
laser  engagement  and  provides  a  means  of  evaluating  the  local  damage. 
Section  IV  develops  the  relationships  necessary  to  convert  the  local 
damage  arising  from  laser  heating  to  damage  parameters  consistent  with 
the  overall  structural  analysis.  Finally,  Section  V  discusses  the 
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results  obtained,  and  Section  VI  presents  conclusions  and  recommenda¬ 
tions  for  further  study. 
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II  Modeling  the  Structure 

The  structural  analysis  method  used  is  based  on  the  displacement 
method  of  finite  element  analysis  (Ref  10,  11,  and  12).  In  such  an 
analysis,  the  continuum  is  replaced  by  a  discrete  model  consisting  of 
a  finite  number  of  elements  connected  at  nodes  (See  Figure  1).  The 
rationale  in  such  an  approximation  is  that  the  response  between  the 
nodes,  i.e.,  the  response  In  the  elements,  can  be  expressed  as  a  func¬ 
tion  of  the  response  at  the  nodes.  Various  interpolation  functions  or 
shape  functions  are  used  to  determine  the  element  response  from  the 
nodal  response.  The  type  of  function  used  depends  on  the  complexity  of 
response  allowed  for  each  element.  The  discretization  reduces  the 
original  differential  equations  of  the  continuum  to  a  set  of  algebraic 
equations  which  can  be  readily  solved  on  digital  computers.  The  itera¬ 
tive  reanalysis  technique  employs  the  original  analysis  of  the  undamaged 
structure.  Therefore  the  development  of  the  reanalysis  technique  will 
also  include  the  development  of  the  basic  steps  of  the  displacement 
method  for  finite  element  analysis. 

The  finite  element  method  used  in  this  study  is  based  on  the  gener¬ 
alized  displacement  method.  The  computer  program  used  is  a  modification 
of  the  program  "ANALYZE"  (Ref  12)  which  was  originally  developed  for  in- 
house  studies  In  structural  analysis  and  optimization  at  the  Flight 
Dynamics  Laboratory,  Wrlght-Patterson  AFB. 

Finite  Element  Method 

The  following  derivation  parallels  a  standard  development  for  the 
finite  element  method  (Ref  10,  11,  12).  The  finite  element  method  was 
derived  through  variational  calculus  using  the  principle  of  minimum 
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potential  energy  as  the  starting  point.  The  potential  energy,  np,  of  a 
structure  is  defined  as  (Ref  11:153) 

np  -  u  +  A  (1) 

where  U  is  the  strain  energy  and  A  Is  the  potential  of  the  applied 
loads.  The  principle  of  minimum  potential  energy  can  be  stated  as 
follows:  "Among  all  displacements  of  an  admissible  form,  those  that 
satisfy  the  equilibrium  conditions  make  the  potential  energy  assume  a 
stationary  value.  Thus 

fill  =  SU  +  6A  =  0 

(2) 

For  stable  equilibrium  Ji  is  a  minimum.  Hence 

P 

S^p  =  fi2U  +  fi2A  >  0  (3) 

Displacements  of  an  admissible  form  are  defined  as  those  satisfying  in¬ 
ternal  compatibility  and  kinematic  boundary  conditions  (Ref  10:56). 

For  a  discrete  finite  element  model  the  total  potential  energy  is  a 
sum  of  the  functionals  for  each  element  (n^),  i.e., 

"<>'  ii  ii <4> 

where  n  is  the  number  of  elements.  The  equations  for  finite  element 
analysis  can  then  be  derived  from  the  element  equations.  If  the  struc¬ 
ture  is  modeled  by  n  finite  elements  connecting  m  nodes,  the  strain 
energy  of  the  element  Is 

l  T  T 
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where  oT  Is  the  transpose  of  the  stress  vector,  is  the  strain  vector 
and  V.|  is  the  volume  of  the  element.  From  Hooke's  Law  the  stress-strain 
relationship  for  a  linearly  elastic  body  can  be  written  as 


?i "  !isi 


(6) 


where  is  the  symmetric  material-stiffness  matrix.  For  the  typical 
homogeneous,  isotropic  plane  stress  element  E^  is  defined  by 


E. 


_ i 

l-v' 


1  0 
v.  1  0 


(7) 


[0  0  \  (l-vi)j 

where  E.  is  Young's  modulus  of  elasticity  and  is  Poisson's  ratio. 

The  finite  element  approximation  is  based  on  the  assumption  that 
the  displacements  within  an  element  can  be  adequately  described  by 
simple  polynomials.  The  coefficients  of  the  polynomials  in  turn  are 
related  to  the  discrete  nodal  displacements  of  the  element.  Therefore, 
the  internal  displacement  equation  is  a  vector  equation  of  the  form 

Mi 


u 


(8) 


where  f^  is  the  vector  of  displacements  in  the  element  coordinate 
system,  is  the  interpolation  or  shape  function  and  u^  are  the 
nodal  displacements. 

The  strain  -  displacement  relations  can  be  written  as 

El  *  §fi  (9) 

where  B  is  a  linear  differential  operator.  For  the  general  problem  B 
is  defined  as 
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Substituting  Equations  6,  8  and  9  Into  Equation  5  gives  the  following 
expression  for  the  element  strain  energy, 

ui '  i  (id 


At  this  point  the  basic  finite  element  assumption,  i.e.,  that  the  in¬ 
ternal  element  displacements  are  a  function  of  the  nodal  displacements, 
makes  its  major  impact  on  the  analysis.  Through  this  assumption  the 
nodal  displacements  have  been  made  Independent  of  the  integration  in 
space  because  they  represent  displacements  at  specific  locations. 
Therefore,  the  can  be  taken  out  of  the  integral  as  follows 

ui  *  r  "I  /  bIiImiM*  s,  Hz) 
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The  elemental  stiffness  matrix,  k^,  is  defined  as 

-kf  '  / 

*1 

Equation  12  can  be  written  as 

U1  ■  7 

The  potential  of  the  applied  loads  is  given  by 


(13) 


(14) 


(15) 


where  is  the  vector  of  lumped  body  forces  per  unit  volume  and  pi  is 
the  vector  of  nodal  forces.  The  sign  is  negative  because  applied  loads 
lose  potential  when  displacement  takes  place.  This  analysis  assumes 
that  the  body  forces  are  zero  allowing  only  nodal  forces.  Thus  Equation 
15  can  be  written  as 


Ai 


(16) 


Substituting  Equations  14  and  16  into  Equation  4  gives  the  following  ex¬ 
pression  for  the  potential  energy 


(17) 


At  this  point,  there  are  n  expressions  which  must  be  summed  to  calculate 
the  potential  energy  of  the  entire  structure.  However,  each  expression 
Is  measured  in  the  local  coordinate  system  of  the  corresponding  element. 
In  order  to  combine  the  element  expressions  into  a  single  structural 
equation,  a  transformation  matrix,  a^,  is  Introduced  such  that 
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“i  *  M 


(18) 


where  D  is  the  structural  displacement  vector  in  the  global  coordinate 
system.  Substituting  Equation  18  into  17  gives 


(19) 


Since  p  is  independent  of  i  it  can  be  pulled  outside  the  summation. 
Then  letting 


and 


(20) 


where  K  is  the  total  structural  stiffness  matrix  and  P  is  the  total 
structural  load  matrix.  Equation  19  becomes 

=  7  PT!Se  -  BTP  (21) 

Taking  the  first  variation  with  respect  to  the  displacements  gives 

snp  =  «(?  PTkd)  ~  «(Btp) 

x  6dTKD  -  SDTP  (22) 

*  «BT(&B  -  p) 
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The  principle  of  minimum  potential  energy.  Equation  2,  requires  that 

55  ■  l  =  0  (23) 

Equation  23  is  a  set  of  algebraic  equations,  which  may  be  solved  numeri¬ 
cally.  Once  the  displacements  are  found  from  the  solution  of  Equation 
23,  the  internal  displacements,  strains  and  stresses  can  be  determined 
by  Equations  18,  8,  9  and  6. 

Finite  Elements 

The  finite  element  program  used  for  this  analysis  contains  four 
basic  elements.  They  are  the  (1)  rod,  (2)  membrane  triangle,  (3)  mem¬ 
brane  quadrilateral  and  (4)  shear  panel.  The  rod  elements  have  constant 
cross-sectional  areas,  and  the  planar  elements  have  constant  thick¬ 
nesses.  The  nodal  displacement  vectors  are  shown  in  Figure  2  for  each 
element.  The  rod  is  a  constant  strain  line  element  allowing  only  axial 
displacements.  The  membrane  triangle  is  a  constant  strain  planar  ele¬ 
ment,  and  the  membrane  quadrilateral  is  constructed  out  of  four  nonover¬ 
lapping  constant  strain  triangles  using  a  fictitious  interior  node  that 
is  removed  using  static  condensation.  The  process  of  static  condensa¬ 
tion  will  be  discussed  in  the  Quadrilateral  Membrane  Subsection.  Con¬ 
struction  of  the  shear  panel  element  is  similar  to  the  membrane  quadri¬ 
lateral  using  four  nonoverlapping  triangles.  However,  only  the  shear 
strains  are  Included  in  the  element  stiffness  matrix  computation.  In 
general,  these  four  elements  are  adequate  for  determining  the  primary 
load  paths  In  most  aircraft  structures  such  as  wings  and  fuselages. 
However,  more  complex  elements  may  be  needed  for  a  detailed  stress 
analysis  of  local  areas. 
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Rod  Element 

The  rod  element  is  an  axial  force  member  commonly  seen  as  a  section 
in  an  analysis  of  a  two  or  three  dimensional  truss  structure.  However, 
within  an  aircraft  structure,  it  is  used  to  model  spar  and  rib  caps. 
Figure  2a  shows  the  allowable  nodal  displacements  as  defined  by  the 
local  coordinate  system.  The  positive  x-axis  is  defined  as  being 
directed  along  the  line  from  the  first  node  to  the  second.  The  rod  is 
a  two  degree  of  freedom  element  since  its  element  displacement  vector 
has  two  components  u^  and  u 2- 

Since  the  rod  is  defined  as  a  constant  strain  element,  the  internal 
displacements  are  assumed  to  be  linear  and  can  be  approximated  by  a 
linear  function 


f  =  c1  +  c2x  =  [1  x] 


=  xc 


(24) 


where  f  ic  the  internal  displacement  at  the  location  x,  and  the  c's  are 
coefficients  to  be  determined.  For  a  linear  elastic  material,  the 
stress  in  the  element  will  also  be  constant. 

Applying  the  boundary  conditions 

f  =  u,  at  x  =  0 

(25) 

f  =  u2  at  x  =  L 

one  obtains 

U1  =  C1  (26) 
u2  *  ci  +  c2L 
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If  one  lets 


u  = 


A  = 


and 


(27) 


c 


then  Equation  26  can  be  written  in  matrix  form  as 


u  =  Ac  (28) 

Solving  for  the  coefficient  vector 

c  =  A'^u  (29) 

and  substituting  into  Equation  24  gives 

f  =  xA_1u  (30) 


Then 


N  *  xA"1 


(31) 
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where  N  is  the  shape  function  describing  the  linear  relationship  between 
the  internal  and  nodal  displacements  (See  Equation  8). 

The  strain  -  displacement  relation  for  axial  strain  can  be  written 


as 

\ 


(32) 


Therefore  B,  as  defined  in  Equation  9,  takes  the  form 


B  = 


(33) 


for  the  axial  rod  element. 


For  the  rod  element  the  stress-strain  relationship  is  simply 

v  -  Ee 


(34) 


Substituting  Equations  30,  33,  and  34  into  Equation  13  gives 


h-/ 

I  »/ 


Vi 


.  Ei 


L-x 


x 

L 


1  -1 
-1  1 


B* °] E< 


/  dv 

V.. 


3x 
0  . 


a 


dV 


(35) 
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Since 


A  A. 


dV 


Vi 


(36) 


where  A  is  the  cross  sectional  area  and  L  Is  the  length  of  the  rod. 
Equation  13  takes  the  form 


».  ■  ^  [!-  1] 


(37) 


for  the  rod  element. 

Triangular  Membrane 

The  basic  planar  element  in  this  program  is  the  membrane  triangle. 
It  is  also  used  in  the  construction  of  the  quadrilateral  membrane  and 
the  shear  panel.  Since  it  is  a  plane  stress  element,  it  can  be  used 
effectively  where  the  primary  loading  is  in-plane  forces,  i.e.,  the  top 
and  bottom  skin  of  aircraft  wings,  flanges  of  I  and  box  beams  subjected 
to  constant  normal  stresses,  and  the  skins  of  sandwich  composite  con¬ 
struction.  However,  this  element  is  not  suitable  when  the  stresses 
through  the  thickness  vary  significantly  (plate  bending).  Inappropriate 
use  of  this  element  will  overestimate  the  stiffness  or  generate  a  matrix 
singularity.  Figure  2b  shows  the  allowable  nodal  displacements  as 
defined  by  the  local  coordinate  system.  The  internal  displacements 
are  assumed  to  be  linear  in  x  and  y  and  can  be  represented  by 
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fx  “  C1  +  C2X  +  c3* 
v  ‘4 +  v  *  c‘y 


where  f  and  f  are  the  x  and  y  displacements  in  the  plane  of  the  plate 
x  y 

measured  in  the  local  coordinate  system  at  location  (x,y)  and  the  c's 
are  the  coefficients  to  be  determined.  Equation  38  can  be  written  in 
matrix  form  as 


f  =  x  c 


(39) 


where 


17 

•  '( 


Since  the  locations  of  the  nodes  are  within  the  domain  of  the  element, 
the  boundary  conditions  are 


}  ■ 


u,  =  c,  +  C2x,  +  c3y, 
.  c4  ♦  c5x,  ♦  Cgy, 

u2  =  C1  +  C2X2  +  c3y2 

v2  =  c4  +  C5X2  +  c6^2 

u3  -  c!  +  c2x3  +  c3y3 

v3  "  c4  +  C5X3  +  c6*3 


(40) 
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where  the  and  are  the  x  and  y  displacements  of  the  1th  node  at 
location  (x.j,  y^)  (See  Figure  2b).  Regrouping  Equations  40  and  re¬ 
writing  in  matrix  notation  gives  the  nodal  displacement  equations  in 
the  form 


Note  that  the  nodal  coordinate  matrix  on  the  right  hand  side  partitions 
into  a  diagonal  matrix.  Since  the  construction  of  the  x  and  y  coeffi¬ 
cients  is  identical,  the  inversion  of  the  partitioned  diagonal  matrix 
is  simply  the  inversion  of  the  component  matrix.  In  matrix  form  Equa¬ 
tion  41  becomes 
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with  the  solution  in  the  form 


For  simplicity,  only  the  derivation  of  x  direction  coefficients  will  be 
shown  and  use  of  the  '  notation  will  represent  the  x  direction  compo¬ 
nent  of  the  partitioned  matrix.  Thus  Equation  41  simplifies  to 


U1 

"1  X,  yf 

V 

u2 

— 

1  x2  y2 

c2 

u3 

x3  >3 

C3 

(42) 


which  has  the  same  form  as  Equation  28  for  the  rod  element.  Solving  for 
the  coefficient  vector 


c  = 


detjAl 


x2y3  -  x3y2 

x3yi  -  xiy3 

xly2  - 

*2  -  *3 

y3  "yl 

yl  - 

*3  '  x2 

X1  ’  x3 

x2  - 

*2 

X1 


(43) 
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where 


Using  the  plane  stress  assumption,  the  differential  operator  B  as 
defined  by  the  strain  displacement  relations  in  Equation  9  reduces  to 


(48) 


which  when  applied  to  N  yields  only  constants.  The  differential  volume 
is  dV^  -  tjdAj  where  t^  is  the  constant  thickness  and  is  the  area 
of  the  element.  E.  is  constant  as  defined  by  Equation  7.  Therefore, 


AFIT/GAE/AA/83M-1 


Equation  13  yields 


-1  *  -Ill-iMiVl 


Substituting  Equations  7,  45  and  48  Into  Equation  49  gives 


V  +x  2^1 

y32  x32  2 


Symmetric 


-(vy32x32  +  x32y32  x32  +  y*2 

Eiti  '(y32y31  +  x32x31  (vX32y31  +  y32x31 

^(1-v?)  (vy32x3i  +  x32y3i  ”^x32x31  +  y32y31 

(y32y21  +  x32x21  ”^vX32y21  +  y32x21 

-(vy32x21  +  x32y21  (x32x21  +  y32y21 


J  +  y2  (1-v) 
y31  X31  '* 

'(vy31X31  +  x31y31  X31  +  y31 

_(y31y21  +  x31x21  ^vX31y21  +  y31x21 

(vy31x21  +  x3]y21  ^j^-)  -(x3iy21  +  y31x2i  ^T"^* 


J  +  x2  LhA 

y 21  X21  2 


'(vy21x2l  +  x21y21  X21  +  y21 


where  Is  the  element  stiffness  matrix  for  a  triangular  membrane 
element  and  x^j  “  xi  “  xj  an<^  y^j  c  y^  -  yj- 
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Quadrilateral  Membrane 

The  quadrilateral  membrane  is  the  element  most  often  used  to  repre¬ 
sent  structural  skins  because  large  areas  with  minimal  curvature  can 
easily  be  represented  by  four-sided  planes.  Figure  2c  shows  the  allow¬ 
able  nodal  displacements  as  defined  by  the  local  coordinate  system.  The 
element  is  assumed  to  be  a  plane  defined  by  the  first  three  nodes.  This 
assumption  ignores  all  warping  and  will  result  in  an  overestimation  of 
the  stiffness  if  out  of  plane  warping  is  significant.  However,  this 
element  can  be  used  in  regions  of  high  warping  if  the  sizes  of  the  ele¬ 
ments  are  reduced  appropriately. 

The  stiffness  of  the  quadrilateral  membrane  is  constructed  by 
dividing  it  into  four  triangular  membranes.  As  shown  in  Figure  3,  a 
fictitious  fifth  node  is  required  and  is  located  by  averaging  the 
coordinates  of  the  element's  four  actual  nodes  using  the  expressions 

X1  +  x2  +  x3  +  x4 
x5  3 

_  y,  +  y2  +  y3  +  *4 

y5  '  4 


This  subdivision  Improves  the  accuracy  of  the  quadrilateral  element  by 
using  more  nodal  displacements  without  burdening  the  user  with  the  task 
of  defining  them.  A  stiffness  matrix,  jjj,  is  computed  for  each  of  the 
four  triangular  elements  using  Equation  50.  Then  the  addition  of  the 
four  matrices  is  accomplished,  similar  to  the  summation  of  the  element 
matrices  described  in  Equations  18  and  20.  A  transformation  matrix,  a^, 
is  Introduced  such  that 
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«1 


*1«Q 


(51) 


where  Uq  is  the  nodal  displacement  vector  for  the  five  node  quadri¬ 
lateral  element,  and  u^  is  the  nodal  displacement  vector  for  the  ith 
triangular  membrane.  Each  triangular  stiffness  matrix  is  then  trans¬ 
formed  and  their  surrmation  results  in  the  following  expression 


1=1 


aT  A  A 


(52) 


where  kg  is  the  10  x  10  stiffness  matrix  for  the  five  node  quadrilateral 
which  includes  terms  for  the  nodal  displacements  of  the  fictitious  node. 
The  fictitious  node's  displacements  are  removed  using  the  following 
static  condensation  procedure  -  a  manipulation  of  the  stiffness  matrix, 
not  an  approximation. 

The  equilibrium  equation  for  the  five  node  quadrilateral  can  be 
written  as 


-kQiiQ  =  Eg 


(53) 


where  Pg  is  the  vector  of  applied  nodal  forces.  Partitioning  Equation 
53  such  that 


*Q 


(54) 
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where  up  are  the  nodal  displacements  for  the  original  four  nodes  to  be 
retained  and  ue  are  those  for  the  fictitious  node  to  be  eliminated  gives 


— 

— 

—  — 

“  — 

Err 

|r 

*re 

Mr 

Br 

>r‘ 

k 

'ee_ 

Be 

_Be 

(55) 


which  can  be  expressed  as  two  separate  equations 

k  u  +  k  u  =  d 
-rr-r  -re-e  cr 


and 


k  u  +  k  u  =  d 
-er-r  -ee-e  "e 


(56) 


(57) 


Since  the  pi  are  the  given  nodal  forces  and  the  fifth  node  Is  not 
actually  present  in  the  original  model,  the  forces  on  this  node  are 
zero.  Therefore,  the  partition  of  p  containing  the  nodal  forces 
applied  to  the  fifth  node,  pg,  must  be  zero.  Using  this  condition  and 
solving  Equation  57  for  ug  gives 


-e  ”  -ee-er^r 


Equation  56  can  now  be  written  as 


(58) 


k_  u 
-rr-r 


*  Br 


or 


-  k  k"^k 
-re*ee-er 


>«r  '  Er 


(59) 
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Note  that  Equation  59  has  the  same  form  as  Equation  23,  i.e..  Equation 
59  is  the  equilibrium  equation  for  the  original  quadrilateral  membrane. 
Thus,  the  element  stiffness  matrix,  k,  can  be  written 


^-rr  ‘  -re-ee-er^ 


(60) 


Shear  Panel 

The  shear  panel  is  a  constant  shear  stress  element  and  cannot  carry 
normal  stresses.  It  was  a  common  practice,  prior  to  the  development  of 
the  finite  element  method,  to  model  wing,  fuselage  and  empennage  struct¬ 
ures  with  constant  shear  panels  surrounded  by  rod  elements  carrying  the 
normal  stresses.  With  the  development  of  finite  element  methods  and  the 
quadrilateral  membrane,  the  top  and  bottom  skins  were  more  easily 
modeled  using  membrane  elements,  since  they  required  the  definition  of 
nodes  and  connectivity  of  a  single  element  rather  than  the  definition  of 
five  elements  needed  for  the  shear  panel  and  rod  construction.  However, 
the  use  of  membranes  to  model  spars  and  webs  results  in  a  gross  over¬ 
estimation  of  the  structural  stiffness  due  to  the  constant  strain 
assumption.  Most  spars,  ribs  and  box  or  I-beam  webs  carry  primarily 
shear  with  some  normal  stress,  resulting  in  deformations  due  to  shear 
stresses,  not  normal  stresses.  In  addition,  the  normal  stresses  usually 
have  steep  stress  gradients  which,  as  stated  in  the  development  of  the 
membrane  triangle,  cannot  be  modeled  using  a  constant  strain  planar  ele¬ 
ment.  Therefore,  the  shear  panel  was  developed  (Ref  13)  to  accurately 
model  these  structural  webs  without  complicating  the  constant  strain 
membranes.  Since  the  shear  panel  is  assumed  to  carry  only  shear 
stresses,  it  must  be  surrounded  by  rod  or  membrane  elements  to  carry  the 
normal  stresses.  Even  though  the  actual  structure  does  not  contain 
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normal  load  members,  the  user  must  provide  fictitious  ones,  usually  rod 
elements,  or  the  stiffness  matrix  will  be  singular.  The  use  of  this 
idealized  model  construction  has  proven  to  produce  accurate  results  when 
used  with  the  appropriate  loading  conditions.  For  example,  when  model¬ 
ing  the  ribs  and  spars  of  a  wing,  the  use  of  shear  panels  is  appropriate 
if  the  major  loading  condition  is  lift  (Ref  12).  However,  when  drag  is 
the  major  loading  condition  shear  panels  will  produce  totally  erroneous 
results. 

The  shear  panel  is  constructed  like  the  quadrilateral  membrane  by 
division  into  four  component  triangles.  However,  the  stiffness  matrices 
of  the  triangles  are  calculated  using  only  the  shear  strain  energy 
terms.  Thus  the  B  matrix  reduces  to 


B  = 


(61) 


which  wh»n  used  in  equation  50  produces 


x 

*32 

Symmetric 

'*32y32 

y32 

Mi 

'*32*31 

y32*31 

x2 

*31 

8A. (l+vi ) 

*32y31 

"y32y31 

'*31y31 

y31 

*32*21 

"y32*21 

'*31*21 

y31 *21 

x2 

*21 

_  *32y21 

y32y21 

*31 y21 

"y31 y21 

'*21y21 

*2 


62) 
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where  =  xi  "  xj  and  =  -  Yj-  The  component  stiffness  matrices 

are  then  combined  and  reduced  in  the  same  way  as  described  for  the  quad¬ 
rilateral  element.  Equations  51  through  60. 

Note  that  since  the  shear  stress  of  an  element  is  dependent  on  the 
orientation  of  the  reference  axis,  the  stiffness  matrix  can  also  be 
affected.  For  rectangular  quadrilaterals,  the  shear  stress  would  be  the 
same  regardless  of  the  side  selected  to  determine  the  reference  axis 
since  all  sides  will  produce  Identical  orientations.  But,  for  the 
general  quadrilateral,  large  differences  in  the  angles  at  the  corners 
can  produce  errors.  However,  by  keeping  the  element  shape  fairly  close 
to  a  rectangle,  errors  will  be  reduced  to  an  insignificant  level. 


Numerical  Solution 

The  majority  of  the  computational  time  for  a  finite  element  analy¬ 
sis  is  spent  solving  Equation  23.  Since  the  matrix  K  can  be  quite 
large,  calculating  its  inverse  could  be  extremely  time  consuming. 
Therefore,  taking  advantage  of  the  fact  that  the  structural  stiffness 
matrix  is  always  symnetric,  positive  definite  and  in  most  cases  sparsely 
populated,  the  solution  can  be  found  economically  using  Gaussian  elimi¬ 
nation  (Ref  14:87-88).  The  procedure  consists  of  three  basic  steps. 

The  first  step  is  a  symnetric  decomposition  of  the  original  stiffness 
matrix  formulated  in  the  following  way 


K  =  L  M  LT 


(63) 
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The  L  and  M  matrices  can  be  calculated  using  the  following  procedure. 
Computing  the  multiplications  in  Equation  63  gives 


Syimietric 

MllLll  +  M22 

M11L31L21  +  M22L32 

M11L31  +  M22L32  +  M33 


MllLnl  i 


MLT- 


M11 

M11L21 

M11L31 


It  follows  then  from  Equation  63  that 


^12  Lln 

«H  =  Kn,  L21  =  . Lnl  ~  RjY 


M22  =  K 


22  "  M11L21 


etc 


The  advantage  of  this  particular  decomposition  is  that  the  sparseness 
characteristic  of  the  stiffness  matrix  is  also  found  in  the  L  matrix. 
This  reduces  the  number  of  computations. 


H  rvi3 
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Substituting  Equation  63  into  Equation  23  gives 


8y  simply  regrouping 


and  letting 


L  M  L'  D  =  P 


L[M  LT  D]  =  E 


Y  *  M  ^  Q 


Equation  66  can  be  written  as 


L  Y  =  P 


The  second  step  is  the  calculation  of  the  vector  Y  by  forward  substitu¬ 
tion.  Rewriting  Equation  67  as 


i  L31  L32 


0  0 

1  0 

Loo  1 


om 


Lnl 


1  M  Pn 


Y1  =  P1 


L21Y1  +  Y2  “  P2 


-  Y2  "  P2  -  l21Vl 


L31Y1  +  L32Y2  +  Y3  =  P3  ^  V  P3  '  L31Y1  "  L32Y2 
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The  last  step  is  to  solve  for  the  nodal  displacement  vector  D  in  Equa¬ 
tion  68.  The  technique  used  is  backward  substitution  which  is  similar 
to  the  forward  substitution  described  in  Equations  70  and  71.  Rewriting 
Equation  68  as 


rw-j-j  o 


0  M, 


22 


1  L21  L31 


0  1 


'32 


0  0  1 


M  0 
nn 


then 


Mnn°n  =  Yn 


Mn-1 ,n-l^n-l  +  Mn-1 ,n-l^n,n-l^n  ~  Yn-1 


n-1 


n-1  M, 


n-1 ,n-l 


'"n.n-l^ 


(73) 


etc 
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Most  of  the  computational  time  for  this  method  is  expended  during 
the  decomposition  (Equation  63),  because  the  forward  and  backward  sub¬ 
stitutions  require  very  little  effort.  Note  that  the  decomposition  is 
independent  of  the  loads  or  the  number  of  separate  loading  conditions. 
Therefore,  an  analysis  using  multiple  loading  conditions  requires  a 
single  decomposition  with  a  forward  and  backward  substitution  repeated 
for  each  loading  case. 

Reanalysis  Techniques 

Repeating  the  entire  finite  element  analysis  for  each  possible  con¬ 
dition  is  so  expensive  it  severely  restricts  the  amount  of  damage  analy¬ 
sis  being  conducted.  Therefore,  the  following  reanalysis  method  has 
been  developed  to  reduce  the  cost,  allowing  comprehensive  vulnerability 
analysis  to  be  performed  (Ref  5).  Assume  the  original  structure  Is 
analyzed  by  the  described  finite  element  program,  then  the  equilibrium 
equation  for  the  structure  is  Equation  23.  When  the  structure  is  dam¬ 
aged,  the  equilibrium  equation  changes  to  the  following  form 

b  5  (74) 

where  K*  is  the  actual  stiffness  matrix  and  p*  is  the  actual  response 
for  the  damaged  case.  The  main  objective  of  a  reanalysis  technique  is 
to  build  an  approximate  solution  of  the  modified  structure,  D*,  using 
the  Information  generated  during  analysis  of  the  original  structure. 
Equation  23  and  the  known  changes.  Therefore,  the  response  of  the 
damaged  structure  can  be  written  in  the  form 

D*  =  D  +  dD  (75) 
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where  the  displacement  vector  D  Is  the  solution  of  the  original  undam¬ 
aged  structure,  and  vector  dD  is  the  true  perturbed  solution  of  the 
damaged  structure.  The  perturbed  solution  can  be  estimated  by  a  trun¬ 
cated  Taylor  series  expansion  as  follows  (Ref  5): 


h  h  a2. 


+  J  JJL—dg.do. 

-  -  -  ii]  39,  ”>  ^,£,^,  39,39,  9' 


where  is  an  implicit  parameter  whose  change  will  affect  the  stiffness 
of  the  ith  element.  There  are  h  ways  in  which  the  stiffness  of  the 
structure  can  be  affected.  The  differential  change  in  response  can  be 


written  as 


h  h  Jlr 


dO  =  l  dg  ♦  L  l  l  -3J—  d9jd9j 

f£l  a9,  '  21  1=1  j=l  s9,39j  9J 


Multiplying  Equation  77  by  the  original  stiffness  matrix  K  gives 


KdD  =  K  T  OL  dg.  +  £.  J  J  dq.dq, 

-  -  '  (£,  39,  S’  Z!  ,£,  j£,  39,39j  °V9J 


Since  the  loads,  P,  are  not  dependent  on  the  g^,  differentiating  the 
original  equilibrium  equation  (Equation  23)  with  respect  to  g^  gives 


*Ld-  -  kIL 

*9i  ‘  ' 


Then  differentiating  Equation  79  with  respect  to  g^  gives 


K  -li..-  «  -  2  1L1L  -  JjL.  0 

’  39i39j  39j  ^  39i 3gj  - 
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Consider  the  last  term  In  Equation  80,  l.e.,  a  f—  in  ’'elation  to 

®i  9j 

the  construction  of  the  stiffness  matrix  K,  defined  by  Equation  20. 
Examine  an  arbitrary  term  from  Equation  20 

$iMi 


Since  the  effects  of  g^  are 
is  independent  of  gj  or 


limited  to  the  jth  element,  the  ith  element 


3ki 


In  addition,  ki  is  assumed  linear  in  gi  giving 


32ki 

.  --J-  =  0 


Therefore,  the  summation  of  Equation  20  gives 


3V9j 


D  =  0 


Substituting  Equations  79,  80  and  83  into  Equation  78  gives 


h  h 


Then  Equation  84  can  be  written  as 


KdD  *  -  dK(D  +  dD) 
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The  iterative  algorithm  for  the  perturbed  solution  can  now  be  written 
from  Equation  85  as 

KdDq+1  =  -  dK(D  +  dDq)  (86) 

where  q  is  the  cycle  of  iteration. 

It  is  very  important  to  note  through  a  comparison  of  Equation  23 
and  Equation  86  that  no  new  analysis  is  required  in  solving  Equation  86. 
Although  the  solution  of  Equation  23  gives  the  unknown  displacement 
vector  D  for  a  known  load  vector  P  while  in  Equation  86  the  objective 
is  to  solve  for  the  perturbed  displacement  vector  dDq+\  the  stiffness 
matrix  K  is  the  same  in  both  equations.  Since  the  stiffness  matrix  was 
already  decomposed  using  Equation  63  for  the  original  undamaged  struc¬ 
ture  it  is  available  for  the  solution  of  Equation  86,  and  does  not  need 
to  be  computed  in  each  iteration.  Therefore  the  solution  of  Equation  86 
reduces  to  the  forward  and  backward  substitution  steps  which  are  the 
least  expensive  of  the  total  solution  steps. 


Estimation  of  the  Effects  of  Yielding 


In  most  everyday  problems  structural  loads  are  designed  to  fall 
within  the  linear  range,  i.e.,  the  loadings  under  consideration 
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generate  stresses  in  each  element  such  that  Young's  modulus  is  a  con¬ 
stant  (See  Figure  4).  However,  with  the  inclusion  of  thermal  effects, 
the  proportional  limit  can  be  reduced  or  the  addition  of  thermal 
stresses  can  put  the  solution  outside  the  linear  regime.  Therefore  the 
reanalysis  program  must  be  able  to  handle  the  effects  of  material  non- 
linearities.  The  most  common  method  for  introducing  material  nonlinear 
ities  into  the  finite  element  analysis  of  structures  is  the  piecewise 
linear  analysis  procedure  (Ref  15).  This  procedure  uses  the  superposi¬ 
tion  of  a  series  of  linear  steps.  The  applied  load  is  divided  into  a 
number  of  small  load  steps,  and  the  stress-strain  curve  is  approximated 
by  a  series  of  straight  lines  (Figure  5)  which  are  used  to  calculate  the 
variation  of  Young's  modulus  with  stress  level.  Each  load  step  is 
applied  sequentially  with  the  Young's  modulus  in  each  step  determined 
from  the  stress  level  in  each  member  calculated  in  the  previous  step. 

The  total  solution  is  obtained  by  adding  the  solution  of  each  step  to 
the  previous  solution,  until  the  full  load  is  reached.  This  procedure 
requires  many  solutions  of  Equation  23  in  the  following  form 

KSADS  =  APS 

where  Ks  is  the  current  stiffness  matrix  using  the  appropriate  Young's 
modulus  for  the  level  of  stress  at  step  s,  aDs  are  the  displacements  for 
step  s,  and  aPs  is  the  load  increment  for  step  s.  Normally  a  minimum  of 
6  to  8  steps  Is  required  for  realistic  results.  For  the  case  of  large 
structures  computing  6  to  8  solutions  is  very  costly.  This  cost  can  be 
reduced  considerably  by  a  modification  of  the  iterative  technique  des¬ 
cribed  in  the  previous  section.  Since  the  nonlinear  effects  consist  of 
Young's  modulus  changes  which  were  previously  Included  in  the  definition 
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STRAIN  (in/in) 


Figure  4.  Typical  Stress-Strain  Curve 
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Fifurc  6.  Stress'-* Strain  Curve  by 
Straight  Line  Approximation 
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of  the  stiffness  change  parameter  gi  in  Equation  82,  the  only  adjustment 
needed  to  the  iterative  technique  is  to  determine  the  amount  of  load  to 
include  with  each  change  of  Young’s  modulus. 

As  before,  the  first  step  in  the  procedure  is  to  make  a  linear 
analysis  of  the  structure  with  the  full  load  using  the  initial  Young's 
modulus.  The  effective  stress  ratios  for  each  element  are  calculated 
using  the  modified  Von  Mises  criterion  (Ref  16)  defined  by 

e{J> 

where  e|^  is  the  effective  stress  ratio  for  the  ith  element  at  the  jth 

break  of  the  stress-strain  curve,  c?x,  and  a ^  are  the  actual  stresses 

in  the  element  and  X.,  Y.  and  Z.  are  the  respective  stresses  correspond- 
J  J  J 

ing  to  the  jth  break  in  the  stress-strain  diagram  (Figure  5). 

For  simplicity  the  following  derivation  will  assume  a  single  break 
in  the  stress-strain  curve.  The  elements  can  then  be  divided  into 
elastic  and  inelastic  groups,  i.e.,  those  whose  stress  levels  fall  with¬ 
in  the  first  and  the  second  straight  line  segments  of  the  stress-strain 
curve,  respectively.  Elements  with  e  <_  1  belong  to  the  elastic  group, 
while  those  with  e  >  1  bilong  to  the  inelastic  group.  Using  the  elec¬ 
tive  stress  ratio,  e,  two  load  factors  are  defined  as 

‘■5s- 

max 

and 

u  =  (1  -  A) 

where  emax  is  the  maximum  effective  stress-ratio  of  all  the  elements. 
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x  is  the  linear  percent  of  the  load  or  that  portion  of  the  stress 
bounded  by  the  break  in  the  stress-strain  curve  and  y  is  the  remaining 
percent.  The  total  applied  load  is  divided  into  two  parts  as  follows 

P  =  XP  +  uP 

The  total  response  can  be  divided  into  corresponding  displacement  vec¬ 
tors  by  simple  scaling  as  given  by 

Dx  =  xD 
Dv  =  pD 

where  Dx  is  the  response  up  to  the  break  and  py  is  the  remaining  re¬ 
sponse.  Since  the  initial  solution  was  calculated  using  the  Young's 
modulus  defined  by  the  first  segment  of  the  stress-strain  curve,  DX  is 
the  valid  solution  for  that  portion  of  the  problem.  However,  Young's 
modulus  changes  In  the  second  portion.  Therefore  the  solution  for  the 
second  segment  can  be  determined  by  the  Iterative  technique  as  follows 

KdDq+1  =  -  dK(Dv  +  dDq)  (87) 

Here  K  is  the  original  linear  elastic  stiffness  matrix  and  dK  is  the 
change  in  the  stiffness  matrix  defined  by 

dS  =  I  n,*, 

where  the  summation  is  over  those  elements  exceeding  the  proportional 
limit,  and  =  (E^  -  Ep/E^  where  Ej  and  E^  are  the  respective  Young’s 
moduli  for  each  segment.  As  described  in  the  previous  section  this 
solution  requires  only  forward  and  backward  substitutions  repeated  for 
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each  iteration  cycle. 

The  total  response  of  the  structure  is  given  by 
D*  =  DX  +  Dy  +  dDq+1  =  D  +  dDq+1 


X  y 


dof’  = 


dcf’ 


where  dDq+^  is  the  displacement  vector  obtained  from  Equation  87,  and 
do<?+^  is  the  stress  vector  in  the  ith  element  due  to  the  displacement 
vector  dDq^ .  In  all  these  computations  it  is  assumed  that  Poisson's 
ratio  does  not  change  during  yielding. 

An  extension  of  this  technique  to  the  case  of  multiple  breaks,  b, 
in  the  stress-strain  curve  is  simply  a  repetition  of  the  procedure  as 
new  breaks  are  encountered.  The  load  factors  for  each  break  are 
computed  by 


(j)  . 


max 


v 


(j)  =  XJ+1 


X 


j 


Using  the  above  definitions  the  solution  for  each  section  can  be  deter¬ 
mined  from  Equation  87.  The  total  response  can  be  written  as 


D*  =  DX  + 


2  2 
(Dy  +  dDy  )  + 


3  3 

(Dy  +  dDv  )  +•••+ 


b  b 
(Dy  +  dDp  ) 
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I 


Convergence  Criterion 

The  solution  converges  If  the  perturbed  solution  can  be  written  as 


dD  =  lim  dDq 
«*»  m 

q  « 


(88) 


where  dD0'  is  obtained  from  Equation  86  (Ref  5).  To  show  the  convergence 
criteria  for  this  solution  technique,  rquation  86  can  be  written  as 

dDq+1  =  -  K-1dK(D  +  dDq)  (89) 

Substituting  the  iterative  results  into  Equation  89  one  obtains 

dDq+"'  =  [I  +  $  +  +•••+  jqJ|D  (90) 

where 


The  superscripts  in  Equation  90  are  exponents  denoting  powers  of  the 
matrix  <j>.  However,  the  exponent  on  the  last  term,  q,  is  equal  to  the 
number  of  iteration  cycles  performed.  From  Equation  90  the  iteration 
converges  only  if  [I  +  <j>  +  +•  •  .+  4>q]  converges.  Note  that  the 

convergence  depends  on  the  nature  of  the  matrix  $.  However,  if  <j>q  -►  0 
as  q  +  »,  the  solution  converges.  This  condition  exists  for  small 
changes  in  the  structural  stiffness  (Ref  9). 

It  is  Important  to  stress  that  assurance  of  convergence  is  not  the 
sole  justification  for  using  this  iterative  technique.  To  be  useful 
this  technique  must  be  considerably  cheaper  than  direct  analysis  of  the 
modified  model  for  the  damaged  structure.  The  computation  time  required 
for  this  method  depends  on  the  number  of  Iteration  cycles  needed  for  the 
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solution,  and  the  rate  at  which  the  solution  converges  depends  on  the 
amount  of  change  in  the  stiffness  matrix  as  shown  in  Equation  89.  A 
measure  of  the  effectiveness  of  this  procedure  could  be  done  by  an 
analysis  of  the  number  of  calculations  performed.  However,  many  vari¬ 
ables  such  as  the  variable  bandwidth  of  the  stiffness  matrix,  the  number 
of  elements  damaged,  the  type  of  elements  damaged  and  the  number  of 
iterations  to  be  performed  are  either  very  problem  dependent  or  very 
hard  to  predict.  Therefore  actual  run  times  were  used  to  e/aluate  the 
effectiveness  of  this  technique.  This  produced  the  estimate  that  the 
run  time  for  approximately  30  cycles  of  iteration  is  equivalent  to  the 
time  required  for  one  direct  analysis  of  the  damaged  structure.  Minor 
damages  require  3-5  iteration  cycles,  medium  damage  requires  5-15 
cycles,  while  major  damage  resulting  in  collapse  of  the  structure 
requires  15-30  cycles  (Ref  5).  These  results  show  that  for  many  damage 
cases  the  iterative  technique  will  be  cost  effective. 

The  criterion  for  convergence  used  in  most  simple  Iteration  schemes 
is  the  amount  of  change  of  the  solution  from  one  iteration  to  the  next. 
However,  In  this  problem  the  solution  must  be  in  a  state  of  equilibrium 
requiring 


where  Aq  is  the  work  of  the  external  forces  for  iteration  q  given  by 

Aq  -  1  PDq 
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and  Uq  is  the  internal  strain  energy  of  cycle  q,  defined  as 


From  Equation  93  the  magnitude  of  Aq+^  -  Aq  depends  on  the  magnitude  of 
the  loads  as  well  as  Dq+^  -  Dq.  Thus  the  rate  of  change  of  the  dis¬ 
placements  is  not  a  valid  criterion  for  convergence.  However,  Equation 
92  could  be  used  as  a  criterion  for  convergence.  Although  the  calcula¬ 
tion  of  Aq  is  inexpensive  requiring  only  the  scalar  product  of  two 
vectors,  the  calculation  of  Uq  is  very  expensive  requiring  the  calcula¬ 
tion  of  the  stresses  and  strains  for  all  the  elements.  Therefore  to 
reduce  the  calculations  necessary  to  check  for  convergence,  the  work  of 
the  external  forces  is  defined  as 

Aq  »  A  +  dAq 

where  A  is  the  work  of  the  external  forces  on  the  original  undamaged 
structure  ard  dAq  =  j  PdDq  is  the  work  of  the  external  forces  resulting 
from  the  perturbed  solution.  Then  the  rate  of  change  in  dAq  can  be  used 
as  an  intermediate  criterion  for  convergence.  When  it  falls  within  a 
predefined  limit,  the  strain  energy  Uq  and  the  work  Aq  can  be  determined. 
If  the  difference  between  these  two  quantities  is  over  a  prescribed 
bound,  the  limit  for  the  rate  of  change  in  dAq  is  reduced,  and  the  iter¬ 
ation  Is  continued.  Otherwise  the  solution  is  considered  to  have 
converged . 
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III.  Modeling  the  Thermal  Effects 

In  order  to  evaluate  the  structural  damage  caused  by  a  laser  strik¬ 
ing  the  structure,  a  method  must  be  developed  to  model  the  thermal 
effects  of  the  beam.  Although  there  are  many  configurations  of  a  laser 
beam  striking  an  aircraft  structure,  which  would  involve  the  incidence 
of  the  beam  and  the  curvature  of  the  structure,  the  scope  of  this  study 
will  be  limited  to  the  basic  situation  of  a  continuous  stationary  laser 
beam  striking  a  flat  plate  at  a  normal  incidence.  The  beam  is  consid¬ 
ered  axially  synmetric  so  that  heat  flux  occurs  only  in  the  axial  and 
radial  directions.  This  study  will  examine  the  effect  that  the  absorbed 
flux  has  on  the  structure.  The  effects  of  the  surrounding  medium  and 
the  absorptivity  of  the  plate  are  taken  into  account  by  properly  adjust¬ 
ing  the  absorbed  flux.  The  effects  of  heat  loss  from  the  heated  plate 
have  also  been  ignored,  since  the  additional  complexity  required  to  in¬ 
clude  realistic  radiation  or  convection  loss  conditions  cannot  be  justi¬ 
fied  by  the  small  impact  they  would  have  on  the  problem  (Ref  17). 

The  thermal  effects  to  be  included  are  the  loss  of  structure  due  to 
melting,  the  induced  thermal  stresses  resulting  from  thermal  expansion, 
and  the  change  in  material  properties  due  to  temperature  dependence. 
Although  other  effects  could  be  included,  these  represent  the  major 
damage  mechanisms  from  a  structural  point  of  view. 

Heat  Conduction  Problem 

The  purpose  of  the  heat  conduction  problem  is  to  determine  the  tem¬ 
perature  distribution  from  which  the  amount  of  melting,  the  thermal 
stress,  and  the  change  In  material  properties  can  be  determined.  Tem¬ 
perature  distributions  are  solutions  to  the  differential  equation  for 
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heat  transfer. 


V 
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where 

Cp  =  specific  heat  (Btu/lbm  °F) 
p  =  density  (Ibm/in  ) 

T  =  temperature  (°F) 

k  =  thermal  conductivity  (Btu/Sec  in  °F) 
t  =  time  (sec) 

Q  =  rate  of  energy  gain  due  to  laser  flux  absorbed 
o 

(Btu/in  sec) 

x,y,z  =  orthogonal  directions 

Consider  the  flat  plate  to  be  a  cylindrical  disk  and  the  prescribed  flux 
to  vary  only  in  the  radial  direction.  Then  no  heat  transfer  occurs  in 
the  circumferential  direction,  and  the  problem  reduces  to  two-dimensions 
in  a  cylindrical  coordinate  system.  Exact  solutions  have  been  found  for 
this  problem  for  constant  properties  and  temperatures  below  melting. 
However,  for  problems  with  a  sufficiently  large  flux  to  cause  tempera¬ 
tures  to  reach  the  melting  point,  available  exact  solutions  are  limited 
to  the  case  of  one  dimensional  problems  subjected  to  a  uniform  flux. 
Therefore,  a  numerical  solution  must  be  used. 

The  solution  could  be  found  using  a  finite  element  analysis  similar 
to  the  technique  described  in  Section  II.  But,  since  the  problem  Is 
time  dependent,  the  solution  would  require  a  complete  analysis  at  each 
time  step.  Similarly,  the  use  of  a  traditional  finite  difference 
approach  also  requires  a  matrix  inversion  at  each  time  increment. 
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Therefore,  the  numerical  approach  used  in  this  study  is  a  modification  of 
an  explicit  technique  that  was  initially  reported  by  Dusinbeere  (Ref  18). 
This  technique  Is  based  on  dividing  the  disk  into  finite  cells  and  per¬ 
forming  a  heat  balance  on  each.  The  advantage  to  this  approach  is  that 
it  requires  less  computer  storage  and  run  time  than  either  of  the 
traditional  approaches  and  yields  comparable  results  (Ref  1). 

Numerical  Solution 

The  first  step  in  this  procedure  is  to  divide  the  disk  into  a  num¬ 
ber  of  finite  cells.  The  thickness,  h,  of  the  disk  is  divided  into  a 
number,  N^.,  of  layers,  and  similarly  the  radius,  r,  is  divided  into  NR, 
sections.  Since  the  structural  model  in  this  study  was  generated  with 
elements  of  a  single  material,  each  cell  is  taken  to  have  the  same  prop¬ 
erties.  For  simplicity,  the  following  development  assumes  that  all 
divisions  are  equal  through  the  thickness  and  that  each  of  the  radial 
divisions  are  of  equal  length,  however,  only  minor  changes  would  be 
required  to  allow  for  unequal  divisions.  Thus  the  disk  Is  modeled  as 
Nc  x  NR  cells  (Figure  6. a).  Each  cell  is  then  referenced  as  a  two 
dimensional  array,  where  the  first  index  refers  to  the  thickness  layer 
and  the  second  refers  to  the  radial  division  (Figure  6.b).  By  using 
this  indexing  scheme  all  properties  associated  with  cell  1,j  can  be 
referenced  in  the  same  manner. 

Consider  a  typical  cell  shown  in  Figure  7.  Cell  i,j  is  surrounded 
by  four  other  finite  cells,  and  heat  conduction  is  the  only  form  of  heat 
transfer.  Using  Fourier's  law  of  heat  conduction 
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where  A  Is  the  area  normal  to  the  direction  of  the  heat  flow  and  x  is 
the  variable  length  between  the  points  of  interest,  the  approximation 
for  the  heat  flow  conducted  out  of  cell  i,j  can  be  given  as 

QRi,j  =  kRi ,jARi ,j 
and 

QZi,j  =  kZi,jAZi,j 

Here  Q  is  the  heat  conducted,  a  the  length  the  heat  travels,  and  the 
subscripts  R  and  Z  denote  the  radial  and  axial  directions  respectively. 

To  perform  a  heat  balance  on  each  element,  it  is  necessary  to  con¬ 
sider  the  location  of  the  element  with  respect  to  the  entire  disk.  If 
an  element  is  located  on  the  surface  of  the  disk,  it  is  possible  that 

energy,  F.  .,  could  be  added  by  the  absorption  of  external  flux  due  to 
■  *  J 

the  laser,  and  there  would  not  be  convected  energy  incoming  from  the 
thickness  direction  since  there  are  no  cells  above  it,  i.e.,  from  a  cell 
of  lower  index  i.  For  the  center  cells,  incoming  conduction  in  the 
radial  direction  is  zero,  since  they  are  solid  disks.  In  addition,  the 
bottom  cells  and  outer  edge  cells  have  no  outgoing  convected  energy  in 
the  thickness  and  radial  directions,  respectively,  due  to  the  imposed 
insulated  boundary  conditions.  These  are  accounted  for  by  setting  to 
zero  the  appropriate  Q.  Considering  each  cell  to  have  flux  values  and 
convection  terms  on  all  four  sides,  and  combining  the  external  flux 
term  with  the  incoming  convected  energy  terms  and  subtracting  the  out¬ 
going  energy  terms,  the  general  expression  for  the  heat  balance  of 


_i>  J  _I _ ^-> 

Ti,j  '  Ti+l,J 
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cell  i,j  is  given  by 


AQ-i 


=  F1,j  +  QRi,j  +  QZi,j  "  QRi ,j  "  QZi,j 


where  aQ.  ^  is  the  amount  of  energy  remaining  to  cause  a  change  'in  the 

1  t  J 

temperature  of  cell  i,j. 

Providing  the  melting  temperature  has  not  been  reached,  the  tem¬ 
perature  rise  in  cell  i,j  during  a  time  increment  At  can  be  written  as 


AQi  At 


where  M.  .  is  the  mass  of  cell  i,j. 
'  »  J 


Melt  Calculations 

In  some  later  time  increment.  At,  the  temperature  will  reach  the 
melting  temperature.  When  the  melting  point  is  reached,  the  cell  has 
not  actually  melted,  and  more  energy  must  be  added  to  bring  about  the 
complete  phase  transformation  of  all  the  mass  in  the  cell.  Thus,  the 
total  energy  required  to  melt  cell  i,j  after  it  has  reached  the  melting 
temperature  is  given  by 


i.j 


r.  .M.  . 

i.j  i,j 


where  r.  .  is  the  heat  of  fusion.  To  model  this  phenomena  an  energy 

'  »  J 

bank  is  established  for  each  cell  with  the  appropriate  y?  4  value. 

Then,  once  the  cell  i,j  reaches  the  melting  temperature,  AtaQ^  is 

subtracted  from  y?  .,  and  the  temperature  of  the  cell  is  held  at  melt- 
*  »  J 

ing.  At  i.  time  increments  after  cell  i,j  reaches  melting  temperature 


k  -J, 
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and  the  entire  cell  is  presumed  to  have  melted  in  the  time  increment 

when  y .  .  -*■  0. 

■  «  J 

Since  the  melted  material  will  add  little  or  no  stiffness  to  the 
structure,  the  melt  is  considered  to  be  instantaneously  removed.  There¬ 
fore,  the  conductivities  of  that  element  are  set  to  zero,  and  any 
external  flux  received  is  transferred  to  the  cell  immediately  below, 
i.e.,  from  cell  i,j  to  cell  i+l,j. 


Laser  Flux  Profile 

The  incident  laser  flux  is  assumed  to  be  circular  with  either  a 


Gaussian  or  uniform  profile.  Neither  of  these  approximations  represents 
the  exact  spatial  profile,  because  experience  has  shown  that  in  reality 
hot  spots  can  occur  (Ref  19:13).  However,  the  Gaussian  profile  is  con¬ 
sidered,  in  general,  to  be  a  suitable  approximation  for  this  type  of 
parametric  study,  because  of  the  random  nature  of  the  intensity,  size, 
and  location  of  such  hot  spots. 

If  the  Gaussian  profile  is  used,  the  incident  energy  absorbed  by  a 
cell  i,j  on  the  surface  is  determined  by 


FoAZi.je' 


.5R2/o 


2) 


where  Fq  is  the  ^  ak  absorbed  intensity,  R  is  the  distance  from  the  cen¬ 
ter  of  the  beam  to  the  center  of  cell  i,j  having  an  area  A-,^  j,  and  2c 
is  the  radius  of  the  laser  beam.  The  assumed  Gaussian  beam  is  actually 
of  infinite  radius,  but  by  defining  the  beam  radius  in  this  manner  86.5* 
of  the  total  energy  absorbed  by  the  plate  is  within  this  radius.  F^  j 
Is  taken  to  be  zero  for  R  -*■  2o,  thus,  the  remaining  13.5%  of  the  energy 
is  neglected. 
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When  the  uniform  profile  is.  selected,  the  incident  energy  absorbed 
is  given  by 

Fi,j  "  FoAZi,j 

for  R  <_  2o  and  F.  .  =  0  for  R  -*■  2o. 

■  •  J 

Results  of  the  Temperature  Distribution 

The  output  of  the  computer  program  is  the  response  of  a  structure 
incurring  a  laser  strike.  Although  the  temperature  distribution  plays 
an  important  part  in  the  final  solution,  it  is  not  possible  to  observe 
only  the  behavior  of  the  temperature  solution  once  it  has  been  integra¬ 
ted  with  the  structures  program.  Therefore  this  section  deals  with  the 
accuracy  of  the  temperature  solutions  and  discusses  the  effects  that 
some  structural  configurations  could  have  on  the  solution. 

To  measure  the  accuracy  of  the  technique,  comparisons  between  the 
temperature  fields  predicted  by  the  numerical  technique  and  exact  solu¬ 
tions  or  other  approximate  solutions  were  made  by  Torvik  (Ref  1).  To 
ensure  the  accuracy  of  the  program  developed  for  this  study,  a  compari¬ 
son  was  made  for  the  one  dimensional  heating  of  an  semi-infinite  solid. 
The  exact  solution,  valid  to  the  onset  of  melting,  is  given  by  (Ref  20) 


T  - 


2Fo(ax) 


1/2 


ierfc 


(88) 


uhere  a  is  the  diffusivity,  k/pCp,  and  ierfc  denotes  the  first  integral 
of  the  complementary  error  function.  The  front  surface  temperature  Is 
given  by 
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T  - 


ax 


1/2 


(89) 


The  comparison  was  made  using  the  properties  of  Al^O^.  This  mater¬ 
ial  was  chosen  as  it  is  a  poorer  conductor  than  metals,  producing  higher 
thermal  gradients  and  a  more  stringent  test  of  the  numerical  process. 

The  properties  are: 

3 

p  =3.8  gm/cm 

k  =  0.104  joule/(cm  sec  °c) 

0^  =  0.885  joule/ (gm  °c) 

Tmel t  =  2313°K 

F  =  4000  joules/(cm  sec) 

Tq  =  300° K 

The  temperature  profiles,  determined  using  Equations  88  and  89,  are 
plotted  at  0.02,  0.04  and  0.06  seconds  in  Figure  8  as  solid  lines.  The 
temperature  distribution  at  the  mid-points  of  each  cell,  determined  by 
the  numerical  procedure  described  in  the  preceeding  section  using  a 
layer  thickness  of  0.015  cm  and  a  time  step  of  0.001  sec,  is  shown  by 
the  asterisks  in  Figure  8.  As  can  be  seen,  the  results  of  the  exact 
and  numerical  solutions  are  indistinguishable. 

In  addition  to  using  the  numerical  solution  to  calculate  known 
solutions  to  verify  the  accuracy  of  the  program,  the  heat  transfer  pro¬ 
gram  was  used  to  observe  the  effects  of  an  internal  structure.  Relating 
the  temperature  distribution  to  an  external  surface  that  can  be  approxi¬ 
mated  by  a  flat  plate  Is  a  fairly  simple  procedure  involving  the  calcu¬ 
lation  of  the  radial  distance  a  given  point  from  the  beam  center  and 
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applying  the  temperature  calculated  at  that  distance.  However,  applica¬ 
tion  to  a  realistic  wing  structure  involves  calculating  the  effects  of 
ribs  and  spars  and  the  temperature  distributions  within  them.  The 
addition  of  structure  below  the  skin  will  act  as  a  heat  sink  and  could 
result  in  a  much  cooler  temperature  distribution  on  the  skin.  To  esti¬ 
mate  the  effect  of  Inner  structure  on  the  temperature  distribution,  a 
"worst  case"  was  examined.  The  "worst  case"  involved  placing  an  in¬ 
ternal  structure  directly  under  the  beam  spot.  The  resulting  tempera¬ 
ture  distribution  was  plotted  against  the  temperature  distribution 
calculated  for  the  flat  plate.  To  observe  possible  effects  on  the 
melting  mechanism,  laser  application  times  of  0.1  and  1.0  sec  were  used, 
giving  one  case  with  minor  surface  melting  and  one  case  with  melt  com¬ 
pletely  through  the  skin. 

For  both  endurance  times  the  inclusion  of  internal  structure  re¬ 
sulted  in  only  slightly  lower  temperatures  within  the  beam  spot.  The 
influence  away  from  the  beam  diminished  with  distance  (See  Figure  9). 

In  addition,  the  area  of  the  plate  removed  by  melting  is  not  affected. 
Therefore  this  study  will  consider  the  effects  of  any  internal  structure 
to  have  no  effect  on  the  temperature  distribution.  This  results  in  a 
conservative  approximation  of  the  damage  incurred,  due  to  the  slightly 
hotter  temperature  solution. 
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Figure  9.  Effect  of  Internal  Structure  on 
the  Temperature  Distribution 
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IV  Modeling  the  Damage 

Once  the  heat  conduction  problem  has  been  solved,  the  damage 
incurred  by  the  structure  must  be  modeled.  The  types  of  damage  to  be 
considered  are:  loss  of  structure  due  to  melt,  addition  of  thermal 
stresses  resulting  from  thermal  expansion,  and  dependence  of  Young's 
modulus  on  temperature. 

Damage  Due  to  Material  Loss 

The  damage  due  to  material  loss  is  defined  as  the  loss  of  stiffness 
that  a  finite  element  incurs  when  some  of  its  material  is  removed 
through  melt.  In  order  to  model  this  phenomenon,  a  parametric  study 
was  conducted  by  modeling  a  flat  plate  in  a  state  of  tension  with  holes 
of  varying  diameters  and  depths.  The  plate  and  the  hole  were  considered 
to  be  symmetric  about  both  the  x  and  y  axis.  Symmetric  boundary  condi¬ 
tions  reduced  the  problem  to  a  quarter  section  of  the  plate  (Figure  10). 
The  symmetric  boundary  conditions  used  restricted  those  points  lying  on 
the  y-axis  from  motion  in  the  x  direction  and  those  points  on  the  x-axis 
from  motion  in  the  y  direction.  The  symmetry  reduction  was  also  consid¬ 
ered  when  the  loads  were  calculated.  To  simulate  a  uniform  load  across 
the  entire  edge,  the  magnitude  of  the  loads  at  each  corner  must  be 
halved,  because  the  mirror  Image  symmetric  section  would  also  have  a 
node  at  that  point  with  a  equal  load. 

When  the  laser  beam  spot  is  considered  to  be  round,  the  material 
loss  will  appear  as  a  circular  hole  partially  or  entirely  through  the 
thickness.  Although  most  of  the  cases  studied  consisted  of  circular 
holes,  non-circular  holes  were  also  considered.  In  an  actual  structure, 
non-normal  Incidence  or  beam  jitter  along  a  preferred  axis  could  produce 
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a  non-circular  beam.  The  resulting  damage  would  then  be  somewhat  oval. 
Therefore,  oval  holes  through  the  entire  thickness  with  the  major  axis 
transverse  to  the  load  were  also  considered  so  that  the  stiffness  loss 
could  be  compared  with  that  for  circular  holes. 

The  stiffness  for  the  damaged  plate  was  calculated  as  a  percent  of 
the  stiffness  for  the  undamaged  plate  by  dividing  the  average  edge  dis¬ 
placement  of  the  complete  plate  by  the  average  edge  displacement  of  the 
damaged  plate.  The  results  are  shown  in  Figure  11,  where  the  cases 
modeling  circular  holes  are  represented  by  solid  lines.  The  cases  cal¬ 
culated  for  oval  holes  through  the  entire  thickness  are  shown  as  aster¬ 
isks.  This  data  reflects  the  effect  of  an  elliptical  versus  a  circular 
shaped  hole  on  the  plate  and  shows  the  stiffness  of  plates  with  oval 
holes  to  be  greater  than  those  with  circular  ones.  Thus  using  the  cir¬ 
cular  hole  values  gives  a  conservative  estimate. 

The  shape  of  the  sides  of  the  hole  must  also  be  considered,  because 
the  hole  is  not  in  reality  a  cylinder  perpendicular  to  the  plane  of  the 
plate.  Such  an  assumption  neglects  a  wedge  of  material  at  the  lower 
edge  of  the  hole.  However,  because  of  the  high  temperature,  the  mater¬ 
ial  in  the  wedge  will  not  add  significantly  to  the  stiffness  of  the 
plate  and  can  be  neglected.  Since  the  program  uses  the  total  energy, 
y-  .,  required  to  melt  cell  i,j  to  determine  the  time  at  which  cell  i.j 

1  »  J 

melts,  the  percentage  of  melting  for  the  cell  can  also  be  determined 
from  Y4  4*  The  program  searches  the  4  looking  for  non-zero  values 

■  *  J  »  t  J 

to  determine  the  amount  of  material  melted  and  the  dimensions  of  the 
hole  created.  The  radial  dimension  of  the  melted  area  Is  taken  to  be 
determined  solely  by  the  top  layer  and  the  depth  of  the  hole  by  the 
center  division.  Each  y^  j  is  tested,  starting  from  the  center  yj  -j  • 
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until  a  value  greater  than  zero  is  found  in  y,  N  .  The  radius  of  the 

l,PIM 

hole  is  then  determined  by  RM  *  Nmrar,  where  nmr  is  the  index  of  the 
first  non-zero  term  radially  and  ar  is  the  length  of  the  radial  divi¬ 
sion.  Similarly  each  j  is  searched  through  the  thickness  giving 
TM  =  N^a^,  where  is  the  index  of  the  first  non-zero  term  through 
the  thickness  and  a^  is  the  thickness  of  each  layer.  The  stiffness  loss 
is  then  determined  from  the  hole's  radial  and  thickness  dimensions  by 
interpolating  from  the  data  for  circular  holes  presented  in  Figure  11. 


Addition  of  Thermal  Loads 

With  the  heating  of  elements  due  to  the  laser  strike,  thermal 
expansion  occurs  and  generates  thermal  stresses.  With  the  addition  of 
the  thermal  load.  Equation  5  becomes 


T  o\ 

2iSiJ 


dV 


(90) 


where  is  the  thermal  strain  produced  by  heating  the  ith  element.  e° 
is  given  by 


where  is  the  vector  of  coefficients  of  thermal  expansion  and  ATi  is 
the  change  in  temperature. 

The  first  term  In  Equation  90  will  follow  the  development  presented 
In  Section  II.  Substituting  Equations  6,  8  and  9  into  the  second  term 
of  Equation  90  gives 

-  /  uTsTbTe,.?«v  (a) 

V1 


65 


.  L 


I 


» 


1 

I 

! 

I 


i 

AF IT/GAE/AA/83M- 1 


Because  all  elements  included  in  this  program  are  isentropic  plane 
stress  elements,  the  thermal  strain  vector  is  given  by 


S? 


a  .AT. 
r  1 


1 

1 

LO. 


(91) 


Multiplying  (a)  by  the  material  property  matrix  gives 


Ei“iATi 

IT 


*iT 


l 

l 

-0. 


(92) 


Since  Equation  92  is  independent  of  volume  and  since  all  elements  in 
this  program  were  constructed  using  a  linear  relationship  between  the 
internal  and  nodal  displacements,  n|bT  is  also  independent  of  volume. 
Therefore  (a)  can  be  written  as 


T.|TpT  Ei“i4Ti 

SiMi  TT^vTT 


or 


If 


T  T  T 

u!n!b!  -4-1  V 

IT  -  vj 


fiV 

n  -  v,) 


(b) 


(93) 
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then  (b)  can  be  written  as 


(c) 


With  this  term  Equation  17  gives  the  following  form  for  the  total 
potential  energy  when  thermal  loads  are  included 


"p  ’  ,!,(?  '  Si*l) 


(94) 


Substituting  the  transformation  matrix  defined  in  Equation  18  and 
letting 


the  total  potential  energy  becomes 

n  =  i  dtkd  -  dtp  -  dt? 

P  2 . 


Taking  the  first  variation  with  respect  to  the  displacement  gives 


«np  =  «DT(KD  -  p  -  t) 


and  when  the  stationary  requirement  is  met  the  equilibrium  equation 
becomes 


KD  -  P  -  j  -  ° 


(95) 
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Combining  the  nodal  forces  and  the  thermal  loads  as 

P‘  ■  -P  *  I  (96) 

and  substituting  into  Equation  95  gives 

KD  -  p*  =  0  (97) 

which  is  exactly  the  same  form  as  Equation  23.  Therefore,  by  using 
Equation  94  to  calculate  the  thermal  loads  and  then  combining  the  two 
load  sets,  all  the  equations  for  the  reanalysis  technique  remain  the 
same  and  can  be  used  with  no  modification.  It  is  important  to  note 
that  the  original  solution  upon  which  the  Iterative  technique  is  based 
must  include  the  affects  of  the  thermal  loads. 

Temperature  Dependence  of  Young's  Modulus 

Although  the  major  loss  of  stiffness  to  an  element  would  be  due  to 
a  loss  of  material  caused  by  melting,  the  effect  on  the  entire  structure 
would  not  be  significant  unless  the  melt  occurred  on  a  major  load  carry¬ 
ing  member.  However,  due  to  the  high  conductivity  of  metals,  a  tempera¬ 
ture  rise  can  occur  over  a  large  section  of  the  structure.  This  can 
reduce  the  magnitude  of  Young's  modulus  in  a  number  of  elements  and 
result  in  a  major  stiffness  change  for  the  entire  structure.  Therefore, 
changes  in  Young's  modulus  can  be  the  major  damage  mechanism  for  the 
structure  as  a  whole. 

To  apply  the  temperature  dependence  affect  to  this  damage  program, 
a  table  of  Young's  modulus  vs  temperature  was  developed.  For  the  sample 
problems  in  this  study  the  Young's  modulus  vs  temperature  curve  for 
T2014  aluminium  from  Mil  Handbook  5,  Figure  12,  was  approximated  in 
tabular  form.  The  temperature  for  each  node  was  calculated  from  the 
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Figure  12.  Young's  Modulus  Dependence  on  Temperature 
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temperature  solution,  and  the  elemental  temrerature  was  determined  as 
an  average  of  the  nodal  temperatures.  Using  that  calculated  tempera¬ 
ture,  an  interpolation  of  the  table  resulted  in  the  new  Young's  modulus, 
El,  for  the  element. 

The  reduction  factor  to  be  used  in  the  iterate ve  technique  is  given 
as 


n 


i 


Ei 
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V  Results 


The  laser  damage  program  developed  in  the  present  study  was  applied 
to  two  different  structures.  These  structures  are  of  two  different 
levels  of  complexity,  one  being  a  simple  two-dimensional  plate  and  the 
other  a  much  more  refined  three-dimensional  wing  structure.  These  prob¬ 
lems  demonstrate  the  ability  of  the  program  to  be  used  for  a  smaller 
localized  analysis  as  well  as  for  a  large  aircraft  structure.  A  com¬ 
plete  description  of  the  two  structures  and  the  reanalysis  results  are 
described  in  this  Section. 

Flat  Plate 

The  first  structure  was  a  flat  plate  clamped  on  one  edge,  subjected 
to  a  uniform  tensile  loading,  and  modeled  with  both  triangular  membranes 
(See  Figure  13)  and  quadrilateral  membranes  ^See  Figure  14).  The  mater¬ 
ial  of  the  structure  was  assumed  to  be  aluminum  with  a  Young’s  modulus 

6  3 

of  E  =  10.0x10  psi  and  a  density  of  p  =  0.1  Ibf/in  .  Six  subcases  were 

defined  by  varying  the  number  and  location  of  elements  exposed  to  the 

laser  strike.  Table  1  lists  the  various  cases  considered  for  this 

structure.  The  laser  size  and  strength  were  varied  to  produce  different 

damage  levels  for  each  case. 

The  overall  size  of  the  plate  was  12.0  in  x  8.0  in  x  .1  in,  and  the 

2  2 

elements  were  uniform  in  size  with  an  area  of  2.0  in  and  4.0  in  for 
the  triangular  membranes  and  the  quadrilateral  membranes,  respectively. 
This  structure  was  selected  because  of  its  simplicity  and  its  size. 

Since  it  is  a  two-dimensional  model,  the  damage  effects  are  easier  to 
understand.  Because  the  plate  is  small  and  the  conductivity  of  aluminum 
is  high,  the  laser  strike  affected  the  temperature  of  all  the  elements. 
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These  attributes  define  a  model  that  Is  easy  to  verify  and  will  show  an 
entire  structure  changing  rather  than  the  effects  of  modifying  a  local¬ 
ized  area.  As  a  result,  the  model  is  a  "worst  case"  example. 

In  order  to  provide  a  baseline  to  compare  the  damage  cases,  the 
solution  to  the  original  undamaged  case  was  used.  The  maximum  displace¬ 
ment  was  0.60  in  for  both  models,  and  the  computer  time  required  was  309 
and  331  10-millisecond  tics  for  the  triangular  membrane  model  and  the 
quadrilateral  membrane  model,  respectively.  The  resulting  deformed 
shapes  are  shown  in  Figures  15  and  16. 

Since  the  temperature  of  all  the  elements  changed,  thermal  loads 
were  generated  for  all  the  nodes.  Some  of  these  loads  were  quite  large. 
For  example,  the  maximum  magnitude  of  the  nodal  loads  generated  from  the 
thermal  effects  in  Case  1  was  13,880.33  Ibf.  With  such  large  loads  the 
panels  can  buckle.  However,  buckling  considerations  are  not  included  in 
this  study. 

The  results  for  the  flat  plate  cases  are  summarized  in  Table  2.  In 
order  to  visualize  the  damage  effects,  the  deformed  shapes  for  these 
cases  are  shown  in  Figures  17  through  22.  The  plots  shown  in  these  fig¬ 
ures  represent  both  reanalysis  and  NASTRAN  results.  The  deformed  shapes 
show  the  expected  responses.  The  figures  show  the  deformed  shape  skewed 
to  one  side  representing  some  in-plane  bending.  Since  the  beam  does  not 
strike  on  the  center  line  of  the  plate,  the  side  of  the  plate  receiving 
the  strike  is  hotter  and  incurs  more  damage.  Thus,  deformation  on  the 
side  of  the  laser  strike  is  greater,  resulting  in  the  skewed  appearance 
of  the  damaged  plate  (Figures  17-22). 

In  the  cases  where  the  damage  appeared  minor,  1.e.»  the  change  in 
the  maximum  displacement  from  the  undamaged  case  was  small,  the 
reanalysis  technique  does  not  seem  to  be  as  efficient  as  in  the  case 
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Figure  15.  Deformed  Shape  Undamaged  Flat 
Plate  —  Triangular  Membranes 
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Table  2  Flat  Plate  Results 
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of  damage  due  to  conventional  weapons  (Ref  5).  However,  since  a  temper¬ 
ature  change  In  every  element  caused  the  appearance  of  one  or  more 
damage  mechanisms  In  every  element,  the  Iteration  extended  over  all  the 
elements.  Thus  the  actual  change  In  the  stiffness  matrix  was  a  major 
reduction,  requiring  an  adjustment  of  every  element.  As  a  result,  the 
stiffness  change  in  each  of  the  reported  cases  is  considered  a  major 
damage  condition,  and  the  reanalysis  technique  proved  to  be  as  effective 
or  better  than  expected.  In  a  larger  structure  the  damage  mechanisms 
would  be  more  localized,  and  the  efficiency  will  become  more  obvious. 
This  Is  shown  with  the  Intermediate  complexity  wing  model  described 
later. 

In  addition  to  evaluating  the  computer  run  times,  the  manhours  to 
set  up  the  models  must  also  be  considered.  The  size  of  the  model  also 
affects  the  manhours.  Altering  the  finite  element  model  to  reflect  a 
material  loss  for  the  flat  plate  structure  is  a  simple  task  requiring 
less  than  10  minutes  clock  time.  However,  to  include  the  Young's  modu¬ 
lus  changes  and  the  thermal  loads  is  much  more  time  consuming.  Since 
NASTRAN  is  a  well  known  finite  element  code,  time  comparisons  will  be 
based  on  the  generation  of  data  required  for  it.  Some  method  must  be 
used  to  generate  a  temperature  distribution  for  each  different  laser 
strike  condition.  The  Young's  modulus  effects  can  be  Incorporated  by 
specifying  that  the  material  is  temperature  dependent  and  including  an 
E  versus  T  table.  However,  in  order  for  the  program  to  use  this  data  or 
to  generate  the  thermal  loads,  temperatures  must  either  be  input  for 
each  node  or  each  element.  Without  the  program  developed  during  this 
study,  the  task  of  converting  the  temperature  distribution  results  to 
node  point  temperatures  would  have  to  be  done  by  hand.  The  time 
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required  to  generate  the  data  In  the  correct  form  Is  about  8  manhours 
for  each  strike.  Note  that  even  if  the  laser  strike  conditions  remain 
the  same,  a  variation  of  the  location  of  the  strike  will  require  this 
conversion  to  be  made.  Therefore,  taking  a  simple  model,  such  as  a  flat 
plate  with  laser  strikes  in  three  locations,  the  program  developed  in 
this  study  shows  a  manhour  savings  of  three  days. 

Intermediate  Complexity  Ming 

The  second  structure  was  a  cantilevered  wing  shown  in  Figure  23. 
This  "Intermediate-Complexity  Wing'1  was  chosen  for  study  as  an  illus¬ 
tration  of  the  application  of  the  program  in  the  preliminary  design  of 
a  lifting  surface.  It  is  a  typical  wing  box  structure  clamped  at  the 
root  and  modeled  using  rods,  triangular  membranes,  quadrilateral  mem¬ 
branes,  and  shear  panels.  The  top  and  bottom  skins  are  modeled  using 
triangular  and  quadrilateral  membranes  and  the  spars  and  ribs  by  shear 
panels  with  rods  providing  the  axial  support.  The  finite  element  model 
has  88  nodes  and  158  elements  (See  Figure  24).  The  applied  loading  con¬ 
dition  was  generated  by  using  simplified  pressure  distributions  repre¬ 
sentative  of  a  subsonic,  forward-center-of-pressure  loading  (Ref  5). 

The  material  Is  assumed  to  be  aluminum  with  the  following  properties: 

E  *  10.5xl06  psi,  v  *  0.3,  p  *  0.1  lbf/in3. 

The  damage  cases  for  this  structure  were  constructed  as  the  analy¬ 
sis  proceeded  by  varying  the  locations  and  the  duration  times  of  the 

laser  strikes  using  a  beam  radius  of  4.0  in  and  a  peak  Intensity  of 
2 

Fq  ■  25  Btu/in  .  The  overall  objective  was  to  observe  the  response  of 
the  structure  as  the  damage  Increased  to  the  level  where  the  structure 
would  collapse.  The  various  cases  are  listed  in  Table  3.  Although  each 
of  these  four  cases  results  in  the  loss  of  skin  from  the  upper  surface. 
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Notes: 

Even  numbered  nodes  are  on 


Figure  24.  Intermediate  Complexity  Wing  Model 
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the  program  can  be  used  to  analyze  laser  strikes  on  any  part  of  the 
structure. 

Results  for  the  wing  cases  are  shown  in  Table  4  and  the  deformed 
shapes  are  shown  In  Figures  25  through  29.  This  structure  is  different 
from  the  flat  plate,  because  it  is  three-dimensional  and  the  lower  skin, 
spars,  and  webs  create  a  pseudo  redundancy  within  the  structure.  In 
addition,  the  overall  size  of  the  structure  is  large  enough,  90  in., 
that  the  temperature  changes  are  more  localized.  Therefore  a  major 
damage  to  an  element,  such  as  complete  removal  through  melt,  and  the 
resulting  temperature  distribution  will  generate  a  localized  damage  con¬ 
dition  and  the  change  in  the  total  structure  will  be  minor.  This 
phenomena  can  be  seen  in  Case  1  where  element  27  has  been  entirely 
removed  by  damage,  and  yet  the  maximum  displacement  for  the  structure 
has  changed  by  only  2.3%  from  the  undamaged  response. 

Case  4  represents  a  collapse  condition.  The  laser  damage  induced 
represents  that  of  a  beam  moving  across  the  wing  at  approximately  center 
span.  The  result  is  a  hole  cut  through  93.75  percent  of  the  chord  at 
that  spanwise  location. 

Figure  30  reflects  the  results  of  all  four  cases.  In  this  figure 
the  total  laser  energy  added  is  plotted  against  the  maximum  displacement 
for  each  case.  The  plot  indicates  that  an  energy  threshold  exists, 
above  which  the  damage  level  increases  significantly.  The  energy  added 
In  Cases  1  and  2  are  below  this  threshold  while  Cases  3  and  4  are  above. 

An  undamaged  analysis  was  performed  to  give  a  run  time  value  to  use 
In  the  operating  cost  evaluation.  The  undamaged  solution  was  computed 
In  3325  10-mil 11 second  tics.  For  Case  1,  with  minor  damage,  the  Itera¬ 
tive  reanalysis  program  required  only  15%  of  the  undamaged  solution  time 
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to  predict  the  solution.  But  for  Case  4,  the  collapse  condition,  the 
Iteration  time  was  360%.  It  must  be  noted  that  because  the  overall  size 
of  the  structure  is  large,  the  number  of  elements  damaged  Is  small  com¬ 
pared  to  the  total  number  of  elements.  In  such  a  case  the  iterative 
procedure  has  shown  to  be  as  effective  as  that  reported  for  damage 
incurred  through  conventional  weapons  (Ref  5). 

When  evaluating  the  manhour  savings,  once  again  the  size  of  the 
structure  is  an  important  consideration.  The  addition  of  a  third  dimen¬ 
sion  further  complicates  the  model  and  increases  the  time  estimation. 

As  with  the  flat  plate,  the  change  in  stiffness  of  the  element  due  to 
melt  can  be  incorporated  into  NASTRAN  data  at  a  relatively  low  cost 
(like  10  minutes).  But  the  temperature  calculation,  by  hand,  would  take 
approximately  3  days.  For  investigating  a  large  number  of  damage  cases 
this  program  could  save  a  considerable  amount  of  time.  For  example,  the 
savings  for  the  four  damage  cases  reported  here  represent  12  days. 
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VI  Conclusions  and  Recownendations 


Based  on  the  results  the  reanalysis  program  with  laser  damage  cal¬ 
culations  Is  capable  of  predicting  the  response  of  a  structure  under 
loads  subjected  to  a  laser  strike.  Due  to  the  numerous  assumptions, 
the  numerical  values  computed  should  be  treated  as  tentative,  but  suited 
to  the  conceptual  design  phase.  From  this  premise,  the  following  con¬ 
clusions  and  reconmendations  are  submitted. 

Conclusions 

1.  The  reanalysis  program  developed  has  met  the  objective  of  this 
thesis  study.  That  is,  it  provided  a  sufficient  method  for  the  pre¬ 
diction  of  structural  Integrity  for  a  structure  encountering  a  laser 
strike. 

2.  The  findings  indicate  that  for  most  minor  and  medium  ranges  of 
damage,  the  reanalysis  technique  is  an  efficient  method  for  analyzing  a 
large  number  of  damage  possibilities. 

3.  When  the  manhours  are  compared  for  the  preparation  of  multiple 
damage  models,  the  reanalysis  technique  provides  substantial  savings 
over  the  generation  of  a  model  for  each  case. 

Reconmendations  for  Further  Development 

There  are  several  directions  that  may  be  taken  at  this  point.  The 
intent  of  this  study  was  to  predict  the  response  of  an  isotropic  air¬ 
craft  structure.  With  the  current  trend  to  use  composites  in  aircraft 
structures,  the  addition  of  the  capability  to  analyze  these  structures 
would  be  a  major  improvement. 
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The  Investigation  of  laser  damage  is  not  limited  to  aircraft  struc 
tures,  but  also  includes  a  significant  amount  of  work  with  missile 
structures.  Including  isoparametric  elements  to  allow  for  the  analysis 
of  this  type  of  curved  surface  is  another  area  of  current  interest. 

Other  improvements  that  are  considered  to  be  possible  modifica¬ 
tions,  but  would  not  have  a  major  impact,  are  including  additional 
structural  considerations  such  as  buckling  and  extending  the  heat 
conduction  solution  to  three  dimensions. 
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A  reanalysis  method  to  analyze  the  strength  of  a  structure  which  has  encoun¬ 
tered  a  laser  strike  is  developed.  The  method  accounts  for  the  following  types 
of  laser  induced  damage;  1)  loss  of  structure  due  to  melting;  2)  change  of  mater¬ 
ial  properties  due  to  temperature  changes;  3)  addition  of  load  due  to  thermal 
stress.  The  program  uses  heat  balance  calculations  over  successive  finite  time 
Increments  on  an  array  of  finite  elements  bisecting  the  laser  beam  spot  to  deter¬ 
mine  the  temperature  distribution.  These  results  are  then  converted  to  struc¬ 
tural  stiffness  parameters  and  a  finite  element  based  reanalysis  method  performs 
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the  structural  analysis.  The  program  was  found  to  give  accurate  results  consis¬ 
tent  with  conducting  a  separate  analysis  for  each  damage  condition  but  with  less 
expenditure  of  computer  time  and  manhours. 
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